Geographic isolation and environmental heterogeneity contribute to genetic differentiation in Cephalotaxus oliveri

Abstract Evaluating the contributions of geographic distance and environmental heterogeneity to the genetic divergence can inform the demographic history and responses to environmental change of natural populations. The isolation‐by‐distance (IBD) reveals that genetic differentiation among populations increases with geographic distance, while the isolation‐by‐environment (IBE) assumes a linear relationship between genetic variation and environmental differences among populations. Here, we sampled and genotyped 330 individuals from 18 natural populations of Cephalotaxus oliveri throughout the species' distribution. Twenty‐eight EST‐SSR markers were applied to analyze population genetics, for the investigation of the driving factors that shaped spatial structure. In addition, we identified the outlier loci under positive selection and tested their association with environmental factors. The results showed a moderate genetic diversity in C. oliveri and high genetic differentiation among populations. Population structure analyses indicated that 18 populations were clustered into two major groups. We observed that the genetic diversity of central populations decreased and the genetic differentiation increased towards the marginal populations. Additionally, the signatures of IBD and IBE were detected in C. oliveri, and IBE provided a better contribution to genetic differentiation. Six outlier loci under positive selection were demonstrated to be closely correlated with environmental variables, among which bio8 was associated with the greatest number of loci. Genetic evidence suggests the consistency of the central‐marginal hypothesis (CMH) for C. oliveri. Furthermore, our results suggest that temperature‐related variables played an important role in shaping genetic differentiation.


| INTRODUC TI ON
The demographic history and population structure can be revealed by detailed characterization of the genetic variation of a species (Lee & Mitchell-Olds, 2011). The genetic variation pattern among populations could result from the complex interaction of natural selection, genetic drift, and gene flow (Eckert et al., 2008;Hamlin & Arnold, 2015). Geographic distance and environmental heterogeneity are two key factors contributing to these processes, and further affect the distribution of allele frequencies and genetic structure in populations (Sork, 2016;Wang, 2013). Under isolationby-distance (IBD), gene exchange among populations is restricted by geographic distance and barriers, resulting in the local accumulation of genetic differences that increase with geographic distance (Bradburd et al., 2013;Wright, 1943). By contrast, under isolationby-environment (IBE), gene flow is limited by strong divergent selection or local adaptation to heterogeneous environment, generating positive correlations between genetic and environmental distances (Wang, 2013). These two patterns are not mutually-exclusive and may co-occur among populations (Manthey & Moyle, 2015). Of course, topography, geographic barriers, and postglacial recolonization may also affect the genetic structure (Cushman & Landguth, 2010;Orsini et al., 2013).
Isolation and fragmentation may influence genetic processes for plant populations within species distribution ranges (Provan & Maggs, 2012). Understanding the processes determining species range limits is important to predict the distributions of species under climate change (Afkhami et al., 2014). However, the determinants of range limits for most species are still poorly understood, mainly because such limits are caused by a complex interplay among a large number of biotic, physical, and historical factors (Roy et al., 2009).
Central-marginal hypothesis (CMH) is one of the primary evolutionary hypotheses to explain the evolution of species' range limits (Sexton et al., 2009). CMH predicts that marginal populations exhibit higher genetic differentiation and lower genetic diversity compared with the central populations (Brown, 1984). When these predictions hold, marginal populations tend to occur in reduced quality habitat resulting in the decrease of effective population size (N e ), which could further reduce gene communication among populations and exacerbate genetic drift within populations. Ecologically marginal populations may be separated not only by distance from the core but also experience different biotic and abiotic environments (Munwes et al., 2010). In addition, CMH is also critical for the management and conservation of marginal populations (Hampe & Petit, 2005). Marginal populations may contain unique adaptive traits to persist under rapid climate change (Sexton et al., 2011). If marginal populations with unique adaptive potential perish in face of human-caused climate change, the species may not be able to effectively expand back to marginal habitats or respond to range shifts. Thus, examining the CMH will confirm whether marginal populations are genetically isolated from central populations, and therefore can inform management decisions for a declining species (Micheletti & Storfer, 2015).
Cephalotaxus oliveri (Cephalotaxaceae) is a dioecious woody shrub or small tree, with the distinctive feature of white stomatal bands on abaxial surface of leaves. As an endemic conifer to southern China, the species has a limited distribution, only appearing in montane regions of Chongqing, southern and western Sichuan, Guizhou, Hunan, north-western Hubei, eastern Jiangxi, northern Guangdong, and south-eastern and north-eastern Yunnan, China.
This plant primarily grows along the sides of streams or valleys at altitudes from 300 to 1800 m and prefers warm and humid habitats (Fu et al., 1999). It contains natural anticancer alkaloid such as harringtonine and homoharringtonine, with high medicinal value (Shi et al., 2010;Xiao et al., 2019). However, its natural population size has significantly declined in recent years due to climate change and overexploitation. Ecological niche modeling of C. oliveri suggested that climate change will cause south-to-north population migration with the contraction of southern regions (Wang et al., 2016).
Several studies led to the development and successful utilization of different markers in conifer species for diversity analyses (Giles-Perez et al., 2022;Wei et al., 2022). However, genetic variation analysis in C. oliveri is still in infancy. ISSRs of C. oliveri have been developed to investigate the genetic variation and identify the correlation between adaptive loci and variation of temperature and precipitation. However, the results of the phylogenetic tree and STRUCTURE analysis were inconsistent. Additionally, the effect of other environmental factors such as soil variables on genetic variation was not considered (Wang et al., 2016). Therefore, research based on a suitable means of detection and more environmental factors is necessary to systematically and comprehensively investigate the population structure and local adaptation of C. oliveri.
To explore the role of geographic isolation and environmental heterogeneity in shaping genetic variation and whether the pattern of genetic variation conforms the prediction of CMH, we used 28 EST-SSR makers to investigate the genetic variation of 330 individuals from 18 natural populations of C. oliveri. The focuses of our work were to (a) assess the genetic variation and population structure of this species, and compare the genetic variation between 12 central and 6 marginal populations; (b) investigate the effects of geographic and environmental factors contributing to genetic variation; (c) identify candidate outlier loci and test for correlation between environmental variables and allelic distribution of outlier loci.

| Plant sampling and DNA extraction
According to the geographic location, isolation, size, and habitat characteristics of C. oliveri populations, 18 populations were collected from the whole distribution region in southern China, including 6 marginal and 12 central populations ( Figure 1; Table S1). The central populations were located in Chongqing, Sichuan, Hubei, and Hunan Provinces, because the result of ecological niche modeling showed that these areas are the highest suitable habitat for C. oliveri (unpublished data). While marginal populations were distributed in Yunnan and Jiangxi Provinces, and isolated from central populations, with small coverage high fragmentation and low densities. Fresh needles were collected from 330 individuals and desiccated with silica gel for subsequent use. Geographic and altitude information were obtained via GPS.
Cetyltrimethylammonium bromide (CTAB) approach modified by Su et al. (2005) was operated for genomic DNA extraction. The concentration and quality of DNA were assessed using microspectrophotometer (ThermoFisher) and 1.0% agarose gel electrophoresis, respectively.

| PCR amplification and genotyping
We selected 28 polymorphic EST-SSRs of C. oliveri previously developed by our laboratory (Liu, Zhang, et al., 2021), and the forward primers of these makers were labeled with two fluorescent dyes (HEX and FAM, Table S2). PCR was performed in 25 μL final solution consisting of 12.5 μL of Taq Mix Kit (Tiangen), 50 ng of genomic DNA, 0.1 μM of each primer, and appropriate ddH 2 O. Amplification procedures were carried out as follows: initial denaturation at 95°C for 5 min, then 30 cycles of 95°C for 30 s, corresponding annealing for 30 s (Table S2), and 72°C for 30 s, with a final extension step at 72°C for 10 min. The alleles with fluorescence labeling were sized using a 3730xl DNA Analyzer (Applied Biosystems: ABI). and analyzed using GeneMapper v.5 (ABI).

| Population genetic diversity
Micro-Checker v.2.2.3 was conducted to check the null alleles (Van Oosterhout et al., 2004). POPGENE v.1.32 was used to assess the genetic parameters of 18 populations (Yeh et al., 1999), including number of alleles (N a ) and effective alleles (N e ), Shannon's information index (I), expected (H e ), and observed (H o ) heterozygosity. The number of private alleles (N p ) was examined using the private_alleles function of R package "poppr" (Kamvar et al., 2014). The inbreeding coefficient (F IS ) and allelic richness (A r ) were calculated using the R package "hierfstat" (Goudet, 2005). The departure from Hardy-Weinberg equilibrium (HWE) of loci in each population was tested by GENEPOP v.4.7.0 (Rousset, 2008), the Markov chain parameters were set at 10,000 dememorization steps, and the numbers of batches and iterations per batch were 20 and 5000, respectively.

| Population genetic differentiation and structure
Analysis of molecular variance (AMOVA) was carried out in Arlequin v.3.5.1.2 to detect genetic variation (Excoffier & Lischer, 2010). The diffCalc function of R package "diveRsity" was run to calculate the pairwise F ST of populations, which then was used to examine the parameter N m according to the formula N m = (1 − F ST ) / 4F ST .
The Bayesian clustering was conducted in STRUCTURE v.2.3.4 to infer the genetic structure of 18 natural populations of C. oliveri (Pritchard et al., 2000). The cluster number (K) ranged from 1 to 18, ten runs were performed for each K, with Markov chain Monte Carlo (MCMC) steps of 100,000 iterations followed by a burn-in period of 100,000 iterations. STRUCTURE HARVESTER was implemented to obtain the optimum K value according to the delta K (∆K) values (Earl & vonHoldt, 2012). Then, CLUMPP and DISTRUCT were conducted to plot the results (Jakobsson & Rosenberg, 2007;Rosenberg, 2004).
GenAlEx v.6.5 was applied to perform principal coordinate analysis (PCoA) (Peakall & Smouse, 2012). Moreover, a UPGMA tree for 18 C. oliveri populations was constructed by using MEGA v. 7.0 with Nei's genetic distance (Hall, 2013).  The soil samples were also collected from 18 populations. The electrical conductivity and pH values of the soil samples were determined as soon as possible after water extraction from fresh soil (soil/water = 1/5 weight/volume) using a potentiometer and a conductivity meter, respectively. The moisture content of soil (fresh or air-dried) was measured by drying at 105°C oven. The airdried soil samples were then smashed, and sieved through 0.2-mm mesh for subsequent analysis. The content of soil organic matter was determined through the potassium dichromate volumetry method, which is divided by 1.724 to obtain organic carbon content. The total nitrogen (N) was evaluated by the Kjeldahl method with a mixture of K 2 SO 4 -CuSO 4 and H 2 SO 4 catalysts. The soil was digested in HNO 3 /HCl/HF (3:1:1) solution, inductively coupled plasma optical emission spectrometer was used to measure the total P, K, Na, Ca, Mg, Si, Al, S, Cu, Zn, Pb, Fe, and Mn. Three soil samples were collected for each population, and each sample was analyzed in triplicate. Finally, 20 soil variables of sampling sites were obtained.
To minimize the multicollinearity effects, we used the vif.cca function in R package "usdm" (Naimi et al., 2014) to check the variance inflation factor (VIF) of the three groups of variables, respectively, and removed the variables with VIF > 10. The final variables set contained 22 environmental variables, including six bioclimatic, seven ecological, and nine soil variables (Table S3).

| Isolation tests
The 22 environmental variables were transformed matrix using PCA. The Euclidean distance between population pairs was calculated using the first five principal components as the environmental distance matrix. The pairwise geographic distance matrix was transformed using the geographic coordinates of populations. The pairwise genetic distance among populations using F ST / (1 − F ST ).
To explore the existence of IBD and IBE, we applied Mantel tests to assess the relationships between genetic distance and geographic or environmental distance. Then, we also used the partial Mantel test to assess the effect of one factor on genetic distance by using another factor as a covariable. The Mantel and partial Mantel tests were implemented using R package "vegan" (Dixon, 2003) based on the Pearson's correlation coefficient method with 10,000 permutations. Finally, multiple matrix regression with randomization analysis (MMRR) was performed to estimate the relative contributions of geography and environment to genetic divergence using the MMRR function with 999 permutations.

| Environmental association analysis
Given that IBD and IBE were considered as useful models for describing the observed genetic differentiation (see Section 3), we performed redundancy analysis (RDA) to identify the relative importance of geography and environment on patterns of genetic variation using rda function of the R package "vegan." To avoid spatial autocorrelation, distance-based Moran's eigenvector maps (dbMEM1-dbMEM10) were applied to represent the geographic variables by transforming the coordinates of each individual using the pcnm function of the R package "vegan". Hellinger transformation was used to transform EST-SSR data into response variables using the decostand function in R. The 22 environmental variables and 10 geographic variables were performed forward selection separately to remove the variables that lacked explanatory power for partitioning using ordiR2step function of R package "vegan". Then, RDA was performed between forward selection variables (explanatory variables) and the EST-SSR matrix (response variables) using rda function. In addition, variation partitioning was also performed using the varpart function in RDA. The significance of all RDA models and constrained axes was detected using the anova.cca function with 999 permutations.

| Identification of outliers and association with environmental variables
We applied two methods to detect loci potentially under selection.
Firstly, a Bayesian approach in BayeScan v.2.1 was performed to detect outlier loci (Foll & Gaggiotti, 2008). The parameters set for BayeScan were the default values. We calculated the q-values to reduce the false discovery rate (FDR), therefore, the locus was recognized as an outlier when q-value was less than .001. In addition, the hierarchical island model (HIM) implemented by Arlequin was used for coalescent simulations, and the heterozygosity and F ST with null distributions were derived, thus estimating the locus-specific p-values. Relationships between outliers (alleles) and environmental factors were analyzed using a linear mixed-effects model (LMM).
The lmer function of the R package "lme4" (Bates et al., 2015) was used to construct LMM, with genetic groups as the random effect while 22 environmental variables (VIF < 10) as the fixed effect. The significance of difference was determined through a likelihood-ratio test using the anova.cca function in R.

| Genetic diversity
The genetic diversity indices for 18 populations of C. oliveri are shown in Table 1 Generally, except for Co164, all other loci significantly deviated from HWE (Table S4).

| Genetic variation and population structure
The STRUCTURE results showed that the ∆K was the highest peak when K = 2, therefore, the optimum number of genetic clusters for populations was suggested as 2 ( Figure

| Genetic comparison of central and marginal populations
We also evaluated the intrapopulation genetic diversity and differ-

| IBD and IBE
The Mantel test detected significant relationships between genetic distance and geographic (r = .670, p = .001) or environmental (r = .391, p = .020) distance, indicating significant patterns of IBD and IBE among populations of C. oliveri ( Figure 4; Table S6). The partial Mantel tests showed that when controlling the influence of another factor, the relationship between genetic variation and geographic distance remained significant (r = .609, p = .001), whereas the association between genetic distance and environmental distance was nonsignificant (r = .183, p = .142; Table S6). In addition, MMRR with the separate effect of geographic or environmental factor on the genetic distance both showed significant correlations, and the contribution of environmental factor had a relatively higher rate than that of the geographic distance (β Geo = .0007, p = .001; β Env = .0743, p = .017; Table S6). Whereas considering geography and environment simultaneously, the result suggested that genetic differentiation of C. oliveri was significantly related to geographic distance (β Geo = .0007, p = .001), while the effect of environmental distance on genetic differentiation was nonsignificant (β Env = .0276, p = .215). and RDA2 were mainly correlated with bio9 and bio13, and bio8 and EVI, respectively (Figure 5b).

| Genetic variation and population structure
The level and distribution of genetic diversity of natural species largely represent their evolutionary potential and adaptive ability to future environmental changes (Gao et al., 2012;Jump & Penuelas, 2005). H e is a primary index to measure gene diversity and has been used to measure the level in several studies (Slatkin & Barton, 1989;Tong et al., 2020). In this study, the average H e and P. massoniana (H e = 0.5483; Mei et al., 2021). Genetic diversity of species may be relative to its characteristics and external factors, such as breeding system, genetic drift, and selection . As a conifer species with scattered distribution, C.

TA B L E 3 Redundancy analyses (RDAs) of Cephalotaxus oliveri.
oliveri primarily grew in montane regions with humid and shady environments. Many isolated populations in fragmented habitats were faced with inbreeding and genetic drift, which was also supported by our results of F IS and HWE analysis. Meanwhile, the long generation time and evolutionary process with a low mutation rate might be responsible for the moderate level of diversity of C. oliveri (Hewitt, 2000). In addition, C. oliveri underwent repetitious climatic changes especially Quaternary glaciations, which would result in loss of distributed range and thus loss of crucial genetic diversity (Hewitt, 2000). Finally, human activities might also lead to decrease in population size.
Based on Wright's (1965) opinion, the genetic differentiation level among populations is high when the F ST > 0.25. The genetic differentiation for C. oliveri was high (0.260), which illustrated that a large proportion of the genetic variation occurred within populations. Our results were consistent with the finding that the genetic variation of perennial cross-pollinated plants mainly existed within the population (Wu et al., 2021;Zhou et al., 2020). In addition, the genetic differentiation of C. oliveri in this study was lower than that determined by ISSR (0.396) and chloroplast sequences (0.639), mainly because of the limited sampling size of populations and differences in molecular markers. The N m of C. oliveri populations was 1.018 (Table S5). As a wind-pollination species, C. oliveri grows in understory and wind only can disperse its pollen over a short distance, making it difficult to genetically communicate with other populations over long distances (Zhang et al., 2005).
Moreover, due to the seed gravity and external factors, most seeds are naturally scattered near the mother plant, which also limits the diffusion distance of C. oliveri seeds. Squirrels and birds can disperse a few seeds, but the dispersal distance is relatively limited (Cain et al., 2000;Lu, 2006).
Results from STRUCTURE, PCoA, and UPGMA indicated two distinct genetic clusters of C. oliveri. The three populations of group 2 were geographically closest and were on the same distribution continuum, while the distance from group 1 was more than 500 km.
The results of IBD showed that the geographic distances between populations had a significant contribution to the genetic divergence of C. oliveri, indicating the fundamental role of spatial distance in the genetic divergence of C. oliveri populations. Wumeng Mountains (WM), about 4200 meters high in elevation, have been illustrated to play a major role in the population divergence of the species (Gong et al., 2011;Zhang et al., 2011). Group 1 and group 2 in this study were divided into different groups on the south or north side of WM, the pair F ST values between two sides of WM were higher than between populations on each side of WM. On the contrary, the gene flow between two sides of WM was lower than between populations on each side of WM. Therefore, the patterns of F ST and gene flow suggested that the WM might act as a physical barrier, interrupting the genetic communication between populations on both sides of the WM, which was also supported by the previous research result of chloroplast DNA .
In addition to geographic distance and barriers, environmental differences may also be the critical factor affecting effective gene exchange. In this study, we also detected the IBE pattern of C. oliveri, indicating the important influence of the environmental difference in the survival and reproduction of individuals. The same signature of the IBE pattern was also found in P. cembra (Toth et al., 2019), P. taeda (Lu et al., 2019), and Pseudotaxus chienii . These results imply that geographic distance and environmental heterogeneity are major factors resulting in the genetic divergence in C. oliveri.

| CMH
The genetic diversity of many species was observed to decline towards the margins of the species range (Lesica & Allendorf, 1995).
Previous review reported that 64.2% of studies (86 of 134) showed the expected decrease in genetic diversity, and 70.2% detected increased differentiation towards range margin (Eckert et al., 2008).
Our results found convincing evidence to support the CMH in C.
oliveri from southern China. In the sampled populations covering its distribution range, we detected a decline in genetic diversity for all genetic indices except two, while the genetic differentiation from the center towards margin range showed an increasing trend. The reasons may include isolation, limited gene flow, increased fragmentation, and genetic drift (Pfeifer et al., 2009). Additionally, a random decrease in the genetic diversity of marginal populations may restrict the evolutionary potential of species, and inhibit the ability of the species adapting to environmental conditions outside the range limits (Blows & Hoffmann, 2005). Similar patterns consistent with the CMH have been detected in Euptelea pleiospermum (Wei et al., 2016), Abies sachalinensis (Kitamura et al., 2020), and Taxus wallichiana , but some species in subtropical China have not found such evidence, such as P. massoniana (Ge et al., 2012), Liriodendron chinense (Yang et al., 2016), and Castanopsis eyrei (Shi et al., 2014). This suggests that different plant species in the same area do not always in accord with genetic patterns of CMH (Wei et al., 2016).
The species distribution model indicated that the marginal habitats for C. oliveri were located in south-western and south-eastern China, while the suitable habitat was at the junction of Hunan and Hubei Provinces (unpublish data), which were predicted as two potential refugia for C. oliveri . Due to low genetic diversity and relative genetic isolation, marginal populations will not expand beyond their current borders, despite the existence of suitable habitats. It is predicted that the distribution range of C. oliveri will shift radically or shrink under global warming, when the majority of marginal populations have low genetic diversity and effective population sizes. Thus, population genetics in combination with ecological niche modeling indicated that marginal populations, especially south-western populations are important for conservation.
Our results may potentially influence management decisions for C.
oliveri in face of climate change.

| Identification of environmental predictors and putatively adaptive loci
Detection of IBD and IBE in C. oliveri by the Mantel test revealed that geography and environment were important for the observed genetic variation. RDA was further implemented to distinguish the separate roles of geographic and environmental factors. The result of RDA suggested that the influence of the environment was more obvious than geography in the genetic differentiation of C.
oliveri. The IBE pattern was also detected in P. chienii  and P. tabuliformis (Xia et al., 2018). Sexton et al. (2014) found that a total of 74% of 70 studies showed the signature of IBE patterns, consisting of 37% only exhibited a significant IBE model but no evidence for IBD. The IBE pattern may arise from natural selection among populations during the local adaptation of these populations to heterogeneous habitat environments (Nosil et al., 2005;Wang & Bradburd, 2014). In this study, differential adaptability to diverse habitats offers a preferable explanation compared with IBD in C. oliveri.
With the 16 environmental variables retained, the soil Fe, Pb, P and K content, and bio8 were found to be important determinants of genetic differentiation through the forward selection of RDA.
Six outliers were identified to be strongly related to environmental variables by LMM. Bio8, altitude, soil pH, soil Mn content, bio12, and bio13 were associated with the most loci. Both the RDA and LMM method showed strong divergent selection signals associated with the temperature-related variable. Temperature is a key factor affecting the phenology and physiology of species, their distribution and range, and structure and dynamics of ecosystems (Nyakatya & McGeoch, 2008). Temperature has been verified to be an important selective driver of adaptive variation in different conifers (Di Pierro et al., 2016;Jia et al., 2020). Moreover, the successful growth and reproduction of plants require various mineral nutrients, which must be absorbed from the soil matrix (DalCorso et al., 2014). Soil Fe, Pb, P, and K content were also important variables to influence the genetic variation of C. oliveri. P and K are essential nutrients for plant growth. Fe is involved in many key cellular metabolic processes, including respiration, photosynthesis, chlorophyll biosynthesis, and synthesis of heme iron and iron-sulfur clusters. However, Pb is a nonessential element and harmful to plants.
In this study, we identified six loci under positive selection and two of which encode TFs. Locus Co224 was significantly associated with soil P, Mg, and Mn content, and its encoded sequence is similar to B-box-type zinc finger protein (BXX proteins 22). Zinc finger proteins are a large family of transcription factors in plants and play a key role in regulating growth, development, and response to abiotic stresses in plants (Khanna et al., 2009). BBX proteins belong to a subfamily of zinc finger proteins that include one or more B-Box domains and are considered to participate in protein-protein interactions (Gangappa & Botto, 2014). BBX genes have been isolated and identified from plants such as Betula luminifera (Dou et al., 2015), Arabidopsis thaliana (Khanna et al., 2009), rice , and Malus domestica (Liu, Dai, et al., 2019). Most of these TFs are related to plant growth, development, and stress resistance. Locus Co274 is located in a region tightly linked to the transcription factor RF2b that has high amino acid sequence similarity with RF2a, and both belong to the bZIP TFs family. The bZIP TFs regulate processes such as light and stress signaling, flower development, and seed maturation (Lindemose et al., 2013;Yang et al., 2019). RF2a and RF2b were confirmed to be predominantly localized to vascular tissues. However, their subcellular localization and the accumulation in different organs were different (Dai et al., 2004). The wheat RF2 subfamily bZIP protein gene TabZIP3 was induced by salt or drought treatment and its overexpression caused improved salt stress tolerance (Guo et al., 2019). Moreover, Co14 is located in one gene of DEAD/DEAH box helicase family, which plays an important role in plant growth and development, disease resistance, and stress tolerance (Li et al., 2008;Nawaz et al., 2018). Several researches confirmed that DEAD-box family was involved in abiotic stress response in plants including rice (Tuteja et al., 2013) and tomato (Zhu et al., 2015). Therefore, these loci may participate in the responses of C. oliveri to different abiotic stresses and play a major regulatory role in the response to stress.

| CON CLUS IONS
This study used EST-SSR markers to investigate the population genetics of C. oliveri, an endemic conifer species in subtropical China.
We found that the existing C. oliveri populations maintain a moderate level of genetic diversity and a high level of genetic differentiation.
Population structure analysis identified two genetic groups. The results also supported the CMH that marginal populations exhibited lower genetic diversity and higher genetic differentiation than central populations in C. oliveri. Furthermore, we detected significant IBD and IBE patterns in this study, and environmental differences had more contributions in shaping population differentiation than geographic distance. We identified six outlier loci that were related to environmental variables, bio8 was associated with the most loci.
Our study gained abundant genetic resources of C. oliveri and provided insights into local adaptation. funding acquisition (lead); methodology (equal); project administration (lead); supervision (lead); writing -review and editing (lead).

ACK N OWLED G M ENTS
We thank Li Liu, Ning Li, Xiaoxian Ruan and Ruonan Wang for advice on samples processing and data analysis.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare that there are no conflicts of interest.